% Monte Carlo Program for the Spatial Probit Model 
% Franzese, Hays, Cook 2014
clear all;
matlabpool open 12;

tic; 
r=100;
tr=110;
t=21; 

%DATA GENERATION

W=[0,0,0,0,0,0,0,0.25,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0;0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0.166666666666667,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0.166666666666667,0.166666666666667,0,0,0,0,0,0,0,0,0;0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0.166666666666667,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0.166666666666667,0,0;0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0;0.333333333333333,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;0.2,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0.2,0,0.2,0,0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0.166666666666667,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0.166666666666667,0,0,0.166666666666667,0,0,0.166666666666667,0,0;0,0,0,0,0,0,0,0,0,0,0,0.2,0.2,0,0.2,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0;0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0.25,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0.166666666666667,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0.166666666666667,0,0,0;0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0.25,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0,0,0.142857142857143,0.142857142857143,0,0,0,0,0,0,0,0,0,0,0.142857142857143,0,0,0,0,0,0,0,0,0,0.142857142857143,0,0,0,0,0,0,0.142857142857143,0,0,0,0.142857142857143,0,0.142857142857143,0,0,0,0;0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0.25,0,0.25,0,0,0,0;0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0.2,0,0,0,0,0,0,0.2,0,0,0,0,0,0.2,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0;0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0.25,0,0,0;0.2,0,0.2,0,0,0,0,0.2,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0;0,0,0.125,0,0,0,0,0,0,0,0.125,0,0.125,0.125,0.125,0,0,0,0,0,0,0,0,0,0.125,0,0,0,0,0,0,0,0,0.125,0,0,0,0,0,0.125,0,0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0.25,0,0;0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0.166666666666667,0.166666666666667,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0.166666666666667,0,0;0,0.2,0,0.2,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0;0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0;0,0.25,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0;0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0.25,0,0,0,0.25,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0.2,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0;0,0,0.166666666666667,0,0.166666666666667,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0;0,0,0,0.25,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0;0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0.166666666666667,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0;0,0,0,0,0,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0.166666666666667,0,0,0.166666666666667,0.166666666666667,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0;0.125,0,0.125,0,0,0,0,0,0.125,0,0,0,0,0,0.125,0,0,0,0,0,0,0.125,0.125,0,0,0,0,0,0,0,0.125,0,0,0,0,0,0,0,0,0,0,0,0,0.125,0,0,0,0,0,0;0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;0,0.2,0,0,0.2,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0;0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0.333333333333333,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0.2,0,0,0,0;0,0,0,0,0,0,0,0,0,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0.2,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0,0,0.25,0,0.25,0,0,0,0,0,0,0.25,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;0,0,0,0,0.166666666666667,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0.166666666666667,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;];


n = length(W);
nt = n*t; 

phi = 0.30; 

i_t = eye(t);
WNT = kron(i_t,W);

i_nt1 = eye(n*(t-1));

top=zeros(n,n*(t-1));
right=zeros(nt,n);

TL=[[top;i_nt1],right];

i_nt = eye(nt); 

rho = 0.25;
rand_mat1 = rand(nt,(r/2));
rand_mat2 = (1 - rand_mat1);
rand_mat = [rand_mat1, rand_mat2];

%Preallocation 
beta_niave=zeros(tr,4);
tstat_niave=zeros(tr,4);
se=zeros(4,tr);
pt=zeros(4,tr);
store_u=zeros(nt,tr); 
store_x=zeros(nt,tr);
store_y=zeros(nt,tr); 
vcov_hat=zeros(4,4,tr);


parfor h=1:tr

XNT = rand(nt,1)*5-2; 

store_x(:,h) = XNT;     
    
seed(:,h) = rng; %seed to retrieve/restore data   

u = randn(nt,1);

store_u(:,h) = u; 

G = (i_nt - rho*WNT - phi*TL); 

ystar = G\(-1.5 + 3*XNT + u);

y=zeros(nt,1);

for i=1:nt
if ystar(i,1) > 0
    y(i,1) = 1;
else
    y(i,1) = 0;
end
end

store_y(:,h)=y; 

z = 1-2*y;
Z = diag(z);

%Niave Probit 
ysl = WNT*y;
ytl = TL*y; 

ysl1 = ysl(n+1:nt,:);
ytl1 = ytl(n+1:nt,:);
XNT_t1 = XNT(n+1:nt,:); 
nt1 = nt - n; 

cons = ones(nt,1); 
X_mat = [cons XNT ytl ysl]; %non-spatial 
[b,dev,stats] = glmfit(X_mat,y,'binomial','link','probit','constant','off');
beta_niave(h,:) = stats.beta;
se_niave(h,:) = stats.se; 


%ESTIMATION

p0 = [-1.5 3 0.01 0.01];
lb = [-5; -5; -0.8; -0.8];
ub = [5; 5; .8; .8]; 
A = [0 0 1 1;
     0 0 -1 -1];
b = [0.9; 0.9];


try

options = optimset;
% Modify options setting
options = optimset(options,'Display' ,'iter');
options = optimset(options,'UseParallel','always'); 
%options = optimset(options,'FunValCheck' ,'on');
options = optimset(options,'PlotFcns' ,{  @optimplotx @optimplotfval @optimplotfirstorderopt });
options = optimset(options,'Diagnostics' ,'on');
options = optimset(options,'HessUpdate' ,'bfgs'); %dfp
%options = optimset(options,'InitialHessType' ,'scaled-identity');
options = optimset(options,'Algorithm' ,'interior-point');
options = optimset(options,'LargeScale' ,'on');
options = optimset(options,'MaxFunEvals' ,10e4);
options = optimset(options,'MaxIter' ,10e2);
options = optimset(options,'TolX', 10e-8);
[p,fval,exitflag,output,lambda,grad,hessian] = ...
fmincon(@(p)sprobit_llf_tscs_SAM_const(p,Z,rand_mat,n,nt,WNT,XNT,TL,i_nt,r),p0,A,b,[],[],lb,ub,[],options);
se(:,h) = sqrt(diag(inv(hessian)));
pt(:,h) = p';
vcov_hat(:,:,h) = inv(hessian); 

catch 
    disp('pos definite error') 
end

trial = h

end
%toc 
matlabpool close;
save('sp_tscs_SAM_phi30_rho25');


true = [-1.5, 3, 0.3, 0.25];
est = mean(pt0'); 
std = std(pt0');
rmse = sqrt((true - est).^2 + std.^2);


output = [est; rmse; std; mean(se_naive); std./mean(se_naive)];

